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ABSTRACT 

The MAGIC collaboration has recently reported the discovery of 7-ray emission from the binary system 
LS I +61° 303 in the TeV energy region. Here we present new observational results on this source in the 
energy range between 300 GeV and 3 TeV. In total 112 hours of data were taken between September and 
December 2006 covering 4 orbital cycles of this object. This large amount of data allowed us to produce an 
integral flux light curve covering for the first time all orbital phases of LS I +61°303. In addition, we also 
obtained a differential energy spectrum for two orbital phase bins covering the phase range 0.5 < <j) < 0.6 and 
0.6 < cj) < 0.7. The photon index in the two phase bins is consistent within the errors with an average index 
T = 2.6±0.2 stat ±0.2 sys . LS I +61°303 was found to be variable at TeV energies on timescales of days. These 
new MAGIC measurements allowed us to search for intra-night variability of the VHE emission; however, 
no evidence for flux variability on timescales down to 30 minutes was found. To test for possible periodic 
structures in the light curve, we apply the formalism developed by Lomb and Scargle to the LS I +61°303 data 
taken in 2005 and 2006. We found the LS I +61°303 data set to be periodic with a period of (26.8±0.2) days 
(with a post-trial chance probability of 10 -7 ), close to the orbital period. 

Subject headings: gamma rays: observations — stars: individual (LS I +61°303) — X-ray: binaries 
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The 7-ray binary system LS I +61°303 is located at a dis- 
tance of ^2 kpc and is composed of a compact object of un- 
known nature (neutron star or black hole) orbiting a Be star in 
a highly e ccentric orbit ( e = .72 ±0.15 ore = 0.55 ±0.05 
following ICasares et al.1 {2005) and iGrundstrom et al.l (120071) 
respectively). 

LS I +61°303 was found to display periodic vari- 
ability in the radio , infrare d, optical, and X-ra y bands 
(Taylo r and Gregorvl IT982L iMarti and Paredesl I 995l 
Mend elson and Mazehl 1 19891 and iParedes et all 1 19971 re 
spectively). 

Th e orbital peri od of the system is 26.4960 days 
long (lGregoryll2002l) . The periastron passage, derived from 
the opti cal spectra, is found to be at phase qb — 0.23 ± 
0. 02 in ICasares et al l J2g 05|) and cf> = 0.301 ± 0.011 
in IGrundstrom et al.l (120071) . adopting a zero-phase at To = 
JD 2443366.775. 

t deceased 
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Radio outbursts are observed every orbital cycle at phases 
varyin g between 0. 45 and 0.95 with a 4.6 years modula- 
tion (Gregory 2002). Radio imaging techniques have shown 
extended, radio-emitting structures with angular extensions 
of ^0.01 to ~0.1 arc-sec, where the radio emission orig- 
inates in a t wo-sided, possib ly precessing, relativistic jet 
{file = 0.6) dMassi et alj|2004l) . These extended radio struc- 
tures have led some authors to adopt the microquasar scenario 
to explain the non-thermal emission in LS I +61°303 (e.g., 
iBosch-Ramon et alj|2006t lBednarekll2006l) . Recent high res- 
olution VLBA measurements show a complex and changing 
morphology different from what is expected for a typi cal mi- 
croquasar je t (se e radio images in iDhawan et al.l e.g., 120061 
lAlbert etal] e.g., (2008b). Furthermore no solid evidence for 
the presence of an accretion disk (i.e. a thermal X-ra y com- 
ponent) has been observed (Chernyako va et ai1l2006l) . This 
seems to favor a scenario in which the non-thermal emission 
in LS I +61°303 is powered by the interaction between a p ul- 
sar and the primary star winds (Marasc hi and Treveslll981l) . 

At higher energies, LS I +61°303 was found to 
be spatially coinci dent with the EGR ET 7-ray source 
3 EG J0241+6103 (iKniffen et all [19971) . Variable emis- 
sion at TeV energies was observed with the MAGIC tele- 
scope (lAlbert et all 120061) an d was recently confirmed by 
VERITAS dAcciari et al.ll2008l) . The system showed the peak 
TeV 7-ray flux at phase <\> ~ 0.65, while no very high-energy 
emission was detected around the periastron passage. 

Here we present new MAGIC telescope observations of 
LS I +61° 303. We briefly discuss the observational technique 
and the data analysis procedure, investigate the very high en- 
ergy (VHE) 7-ray spectrum during the high emission phase of 
the source, and put the results into perspective with previous 
VHE 7-ray observations of this system. Finally we analyzed 
the temporal characteristics of the TeV emission and find a 
periodic modulation of the signal with the orbital period. 

2. OBSERVATIONS 

The observations were performed from September to De- 
cember 2006 using the MAGIC telescope on the Canary is- 
land of La Palma (28.75°N, 17.86°W, 2225 m a.s.l.), from 
where LS I +61°303 is observable at zenith distances above 
32°. The telescope operates in the energy band from 50 — 
60 GeV (trigger threshold at zenith angles less than 30 de- 
grees) up to tens of TeV, with a typical energy resolution of 
20 — 30%. The accuracy in reconstructing the direction of 
the incoming 7-rays is about 0.1°, depending on the energy. 
A detail ed description of th e telescope performance can be 
found in lAlbertetaD(l2008cl) . 

The data on LS I +61°303 were taken between 15 th of 
September 2006 and 28 th of December 2006 covering 4 or- 
bital periods of the system. In total 120 hours of data were 
taken at zenith angles between 32 and 55°, with ~ 97% of 
the data below 44°. After pre-selection of good quality data 
a total of 1 12 hours of data remained for the analysis. About 
17% of these were recorded under moderate moonlight con- 
ditions. Due to the different observation conditions such as 
bad weather, too bright moon or too large zenith angle, the 
data set was not uniform with the orbital phase. In Table Q]the 
observation times of the analyzed data are summarized. 

The observations we re carried out in wobble 
mode dFomin et al.1 1 19941) . i.e. by alternately tracking 
two positions at 0.4° offset from the actual source position. 
This observation mode allows for a reliable background 
estimate for point like objects such as LS I +61°303. 



3. DATA ANALYSIS 

The data analysis was carried out using th e standard 
MAGI C analysis and reconstruction software ( Albe rt et al.l 
2008c) . The images were cleaned by requiring a minimum 
number of 10 photoelectrons (cor e pixel s ') and 5 photoelec- 
trons (boundary pixels), see e.g. iFeganl d 1997b . The qual- 
ity of the data was checked and bad data such as accidental 
noise triggers or data taken during adverse conditions (very 
low atmospheric transmission, car light flashes etc.) were re- 
jected. Fr om the rema ining events, image parameters were 
calculated (lHillaslll985l) . 

For the 7/hadron separation a multidimensional clas- 
sification ^grocedure based on the Rand om Forest 
method (lAlbert et al.l l2008at iBock et all 120041) was used. 
For every event a parameter called hadronness (h) is derived, 
based on the values of the events image parameters. The 
hadronness denotes the probability that an event is a hadronic 
induced (background) event. The final separation was 
achieved by a cut in h which was determined by requiring 
80% of the simulated Monte Carlo (MC) 7-ray events to 
be kept. In addition to the cut in h a. geometrical cut in the 
squared angular distance of the assumed source position to 
the shower direction axis (6 2 cut) was performed so that 70% 
of all simulated MC 7-ray events from a point-like source 
are left after the cut. The cut efficiencies were determined 
by optimizing the significance of a Crab Nebula data sample 
recorded under the same zenith angle as the LS I +61°303 
data set. The same cut procedure was applied to the final 
LS I +61°303 sample. The energy of the primary 7-ray 
was reconstructed from the image parameters using also a 
Random Forest method leading to an assigned estimated 
energy for each reconstructed 7-ray event. The differential 
energy spectrum is unfolded t aking into accou nt the full 
instrumental energy resolution (lAlbert et al.l 120071) . For the 
integral flux calculation of the light curves we used fixed cuts 
(for all energies) in hadronness and 9 2 . In the case of the 
energy spectrum determination we derived fixed hadronness 
and 9 2 cuts for each energy bin. 

The main contributions to the systematic error of our anal- 
ysis are the uncertainties in the atmospheric transmission, the 
reflectivity (including stray-light losses) of the mirror and the 
light catchers, the photon to photoelectron conversion cali- 
bration and the photoel ectron collection effi ciency in the pho- 
tomultiplier front-end (lAlbert et al.ll2008ct) . Also MC uncer- 
tainties in the detector simulation and systematic uncertainties 
from the analysis methods contribute significantly to the over- 
all error. 

All errors in this paper are statistical errors, otherwise it is 
stated explicitly. In addition there is a 30% systematic uncer- 
tainty on flux levels and 0.2 on the spectral photon index. 

MAGIC has the capability to operate under moderate 
moonlight. This permits to increase the duty cycle by up to 
28%, thus considerably improving the sampling of transient 
sources. In particular, 17% of the data used in this analy- 
sis were recorded under moonlight. The nights which were 
partly taken under moonlight conditions are labeled with a 
star in Tabl^U For these days we estimate an increased sys- 
tematic error of ~ 40% instead of the ~ 30% in the case of 
the dark night observations. All spectra are derived from data 
which is only taken under dark night conditions and thus no 
additional error is present in the obtained parameters. 

3.1. Light Curve 



Periodic VHE 7-ray emission from LS I +61°303 
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FIG. 1.— VHE (E > 400GeV) gamma-ray flux of LS I +61°303 as 
a function of the orbital phase for the four observed orbital cycles (4 upper 
panels) and averaged for the entire observ ation time (lowerm ost panel). In the 
lower most panel the previous published ( Albert et al. 2006) averaged fluxes 
per phasebin are shown in red too. Vertical error bars include la statistical 
error. 

Figure[T]presents the gamma-ray flux above 400 GeV mea- 
sured from the direction of LS I +61° 303 as a function of the 
orbital phase of the system for the 4 observed orbital cycles. 
The probability for the distribution of measured fluxes to be a 
statistical fluctuation of a constant flux (obtained from a x 2 fit 
to the entire data sample) is 4.4 x 10~ 6 (x 2 /dof = 108.9/51). 
In all orbital cycles significant detections (S > 2a) occurred 
during the orbital phase bin 0.6-0.7. The highest measured 
fluxes are dominantly found in this phase bin. Among those 
nights around phase 0.65, the night MJD 54035.1 1 shows the 
maximum flux, with statistical significance of 4.5 a. 

At the peria s tron passage (phase 0.23, according 
to ICasares et al.l l2005h the flux level is always below 
the MAGIC sensitivity and we derive an upper limit 
with 95% confidence level of 4 x 10~ 12 cm _2 s _1 (MJD 
53997). If we take for th e periastron passa g e the phase 
value 0.3 as obtained by iGrundstrom et al.l (120071) . we 
detect a marginal signal on MJD 53999 with a flux of 
F(E > 400 GeV) = 5.3 ± 2.4 x lO-^cm^s- 1 . Since the 
correct value for the periastron passage is yet debated we can 
put strong constrains to the emission only in the case of phase 
0.23. 

Summing up all data between phase 0.6 and 0.7, where the 
maximum flux level is observed, we determine an integral flux 
above 400 GeV of 

F(E > 400GeV) = (7.9±0.9 stat ±2.4 syst ) x lO^cm-V 1 . 
The above quoted flux corresponds to 7% of the integral 



TABLE 1 

Observation time, orbital phase, integral flux (above 
400 GeV), flux upper limit at the 95% c onfidence level (g iven 
in case flux significance is < 2a, drolke et al.i20051 
following)). Nights partly taken under moonlight 
conditions are labeled with a star. 



Middle Time Obs. Time Phase 



Flux 



Upper limit 



(MJD) 


(min) 




io- 12 

(cm -2 s - 1 ) 


io- 12 

(cm -2 s - 1 ) 


53993.18* 


137 


0.08 


3.7 ± 2.3 


8.5 


53994.17* 


1 12 


0.11 


0.6 ± 2.7 


6.2 


53995.17* 


157 


0.15 


-2.0 ± 2.2 


3.0 


53997.15 


229 


0.23 


0.3 ± 1.8 


4.0 


53998.15 


21 1 


0.26 


2.0 ± 2.0 


6.0 


53999.10 


133 


0.30 


5.3 ± 2.4 




54001.12 


82 


0.38 


-3.6 ± 3.8 


5.1 


54002.09 


188 


0.41 


2.4 ± 2.3 


7.1 


54003.08 


144 


0.45 


1.8 ± 2.7 


7.2 


54004.08 


158 


0.49 


-4.0 ± 2.5 


2.5 


54005.07 


155 


0.52 


3.0 ±2.5 


8.1 


54006.07 


162 


0.56 


1.8 ±2.7 


7.2 


54007.08 


139 


0.60 


4.4 ± 2.8 


10.2 


54008.07 


152 


0.64 


8.8 ± 3.1 




54009.08 


147 


0.68 


4.4 ± 2.6 


9.7 


54013.24 


7 


0.83 


0.8 ± 10.7 


26.7 


54022.10* 


186 


0.17 


1.7 ± 2.0 


5.8 


54023.10* 


269 


0.20 


-2.9 ± 1.5 


1.4 


54024.08* 


20 


0.24 


-0.4 ± 7.0 


15.2 


54029.02 


134 


0.43 


-1.1 ± 2.5 


4.1 


54030.01 


161 


0.47 


-0.4 ± 2.3 


4.2 


54031.01 


163 


0.50 


5.9 ± 2.6 




54032.01 


139 


0.54 


3.4 ± 2.9 


9.2 


54035.11 


150 


0.66 


12.7 ± 2.9 




54039.09 


93 


0.81 


-1.4 ± 1.2 


1.7 


54055.97 


181 


0.45 


4.0 ± 2.2 


8.5 


54056.96 


223 


0.48 


-0.2 ± 2.1 


4.2 


54057.90 


66 


0.52 


3.3 ± 3.8 


11.2 


54058.90 


57 


0.56 


2.3 ± 3.3 


9.2 


54060.00 


17 


0.60 


16.5 ± 6.8 




54061.96 


221 


0.67 


5.9 ± 2.2 




54062.96 


228 


0.71 


5.5 ± 2.1 




54063.95 


56 


0.75 


3.6 ± 4.0 


12.1 


54065.00* 


71 


0.79 


4.5 ± 3.8 


12.4 


54066.02* 


185 


0.82 


1.1 ± 2.4 


5.9 


54067.04* 


188 


0.86 


0.3 ± 2.3 


5.0 


54068.08* 


77 


0.90 


-1.5 ± 3.6 


6.1 


54081.89 


17 


0.42 


-0.3 ± 5.4 


12.6 


54082.85 


77 


0.46 


2.9 ± 3.7 


10.6 


54083.88 


31 


0.50 


4.4 ± 5.3 


15.9 


54084.85 


63 


0.54 


1.5 ±4.5 


10.8 


54085.95 


111 


0.58 


2.4 ± 1.4 


5.5 


54086.95 


282 


0.61 


8.6 ± 1.8 




54088.01 


82 


0.65 


9.7 ± 3.6 




54088.95 


83 


0.69 


3.4 ± 2.9 


9.4 


54089.89 


29 


0.73 


0.4 ± 3.7 


9.0 


54090.88 


176 


0.76 


3.6 ±2.2 


8.1 


54091.90 


140 


0.80 


1.9 ±2.8 


7.6 


54092.92 


92 


0.84 


15.6 ± 3.8 




54093.97 


92 


0.88 


7.0 ±3.5 




54095.01* 


57 


0.92 


1.1 ± 1.1 


4.1 


54096.02* 


49 


0.96 


3.6 ±4.4 


12.8 



Crab nebula flux in the same energy range. The mean flux 
for all other phase bins can be found in Table [2] This is well 
in agreeme nt with the flux m easured by MAGIC in the first 
campaign (Alb ert et al.ll2006l) . The data we presented here 
have been reanalized with an improved energy estimation. 
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TABLE 2 

Average flux level above 400 GeV for each orbital phase bin. 

Flux upper limit at the 95% confidence level are quoted in 
case flux significance is < 2(7 (following rolke e t al.i | |200"5|) 1 . 



TABLE 3 
SPECTRAL FITTING PARAMETERS 



Phase bin 


Flux 


Upper Limit 




(10- 12 cm- 2 s" 1 ) 


(10- 12 cm" 2 s" 1 ) 


0.0-0.1 


3.7 ± 2.3 


8.5 


0.1-0.2 


0.2 ± 1.2 


2.7 


0.2-0.3 


0.3 ± 0.9 


2.2 


0.3-0.4 


-1.2 ±2.8 


4.3 


0.4-0.5 


0.7 ± 0.8 


2.4 


0.5-0.6 


3.1 ± 1.0 




0.6-0.7 


7.9 ± 0.9 




0.7-0.8 


4.3 ± 1.2 




0.8-0.9 


2.8 ± 1.1 




0.9-1.0 


0.7 ± 2.0 


4.8 



One very interesting peculiarity found in the light curve is 
that a second peak in the flux level is seen in the last observed 
period at phase 0.84. The flux is (16 ± 4) x 1CT 12 cm^s" 1 
which is at a similar level compared to the maximum flux 
detected in phase -0.65 (MJD 54035.11). This high flux 
was not seen in similar phases in any previous cycle, where 
only upper limits could be set (see Table Q]). Six hours 
after our measurement, the data of the Swift X-ray satel- 
lite s howed a high flux (0.25 ± 0.01 counts s _1 )at phase 
0.85 dEsposito et al.l l2007). These Swift observations did not 
cover the same orbital phase in any other orbit. Beside this 
second peak the main X-ray emission p eak is found between 
the phases 0.5-0.8 (~ 0.24 counts s -1 lEsposito et"aT1l2007l) 
in exactly the same phase bin where LS I +61°303 is detected 
by MAGIC. This is a hint for a correlated X-ray/TeV emis- 
sion. 

Our measurement i s in agreement with the published VER- 
ITAS measurements (lAcciari et al.ll2008l) . that LS I +61° 303 
is detected at TeV Energies in the phase range 0.5-0.8. The 
MAGIC and VERITAS data are not strictly simultaneous 
taken and VERITAS did not observe LS I +61°303 in De- 
cember 2006 were the second peak occurred. 
Due to our long observation time and dense sampling of or- 
bital phases we obtained the currently most detailed light 
curveofLSI+61°303. 

3.2. Spectral studies 

As seen from the light curve (Fig. [0 LS I +61°303 is a very 
variable source which shows high flux levels only at some 
orbital phases. For the phases 0.5 < <j> < 0.6 and 0.6 < cf> < 
0.7, where we measured significantly high flux levels we were 
able to determine differential energy spectra. In both cases the 
obtained energy spectra are compatible with pure power laws. 
In the case of the phase bin 0.6 < cj) < 0.7 a power law fit 
gives: 

-2.6±U.2 slal ±0.2„ 



dF 
dE 



(2.6 ± 0.3 stat ± 0.8 



systy 



10 



-12 



TeV cm 2 s 



E 



ITeV 



with a reduced y^/dof = 5.22/5. The spectral fit pa- 
rameters a gree excellent w ith the previous reported ones by 
MAGIC (lAlbertetalJl2"006h . 

In addition we derived the differential energy spectra for 
the two nights with a signal of > 4.5er significance, which 
are part of the same phase bin 0.6 < cj) < 0.7. Both spectra 
are also well described by a pure power law (see table[3]). No 
evidence for spectral variations has been found. 



Phase / MJD 


Flux 


Spectral Photon Index 




(lO-^TeV-^m" 2 s" 1 ) 




0.6-0.7 


2.6 ± 0.3 


2.6 ± 0.2 


0.5 - 0.6 


1.2 ± 0.4 


2.7 ±0.4 


54035 


3.6 ± 1.1 


2.7 ±0.5 


54086 


3.2 ± 0.6 


2.6 ± 0.3 




FIG. 2. — Shown is the spectrum of phase bin 0.5 < (j> < 0.6 (green), 
phase bin 0.6 < <f> < 0.7 (red) and the spectrum obta ined from previous 
MAGIC measurements (blue dashed) Albert et al. 12006). The spectra from 
phase bin 0.6 < <f> < 0.7 and the previous MAGIC measurements can be 
well described by a simple power law with photon index 2.6 ± 0.2. The 
spectral slope of phase bin 0.5 < 4> < 0.6 is compatible with these results 
within the errors. 

In case of phase bin 0.5 < (j> < 0.6 we obtained: 

dF = (1.2±0.4 stat ±0.3 syst )-10- 12 / E y 2 - 7±0A ^ 
dE TeV cm 2 s \1 TeV J 

with a reduced x 2 /dof = 1.42/4, showing that the spectral 
shape is well compatible with a simple power law. 

The energy spectra of the two phase bins together with the 
power law fits are shown in Fig. [2] The corresponding fit pa- 
rameters and their errors are also shown in Table [3] 

The spectral indices of the fitted power laws, for both phase 
bins and the single night spectra, are compatible within their 
errors indicating that no significant spectral changes happened 
between the different phase bins and between the different 
orbital cycles. In the phase bins 0.0 < (j> < 0.5 and 0.7 < (j> < 
1.0 the 7-ray flux is too low to derive meaningful differential 
energy spectra. 

Another possibility to search for spectral variation is by 
means of the hardness ratio HR which we define as the ratio 
of the integral flux between 400 GeV and 900 GeV and above 
900 GeV. The HR plotted against the total integral flux above 
400 GeV for each night with a signal above 2cr significance 
is shown in Fig. [3] The requirement of the 2a significance 
of the signal is to minimize systematic effects on the calcu- 
lation of the correlation coefficient. We do not find any clear 
correlation between the HR and the flux level. Thus we do 
not find any change in the spectral behavior in nights where 
LS I +61°303 is detected at modest significance. 

The spectral studies on LS I +61°303 exhibits that the 
spectrum is soft (compared to other galactic sources) dur- 
ing all phases in which LS I +61°303 is detected at TeV 
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FIG. 3.— The HR, defined as F(E > 900GeV)/F(400 < E < 
900 GeV), vs. the integral flux F(E > 400 GeV). Each point is one night 
with a reconstructed signal of at least 2cr significance. There is no clear cor- 
relation between the HR and the flux level. 
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FIG. 4. — Intra-night integral flux behavior for the longest observed night 
18 th December 2006, MJD 54086.95 and orbital phase 0.61. The light curve 
is fitted to a constant flux level with a probability of 90% (x 2 /dof = 1 . 08/4). 

energies. While the flux level changes on timescales of 
days and reaches a maximum flux (detected above 3<r) of 
15.6 x 10~ 12 cm~ 2 s _1 the source shows a spectral photon 
index of 2.6 ± 0.2, compatible with being constant. 

4. TIMING ANALYSIS 
4. 1 . Search for intra-night variability 

LS I +61°303 was found to be variable on timescales 
of days in the previou s observational campaign by 
MAGIC (lAlbert et al.ll2006l) . A still open question is whether 
LS I +61°303 also shows variability on shorter time scales. 
We investigated the data for all nights with a significant flux 
level (F(E > 400GeV) > 4 x lCT 12 cm^s" 1 ) with respect 
to intra-night variability on time scales ranging from 30 to 75 
minutes, in steps of 15 minutes. This yields 16 suitable nights. 
Among those, the longest one was MJD 54086.95 (phase 
0.61). Its intranight light curve is shown in Fig. [4] where each 
bin has ~ 1 hour width. The light curve is fitted with a con- 
stant flux level with probability of 90% (x 2 /dof = 1.08/4). 
For the rest of the nights and tested time scales, the post-trial 
probabilities of being chance fluctuations of a constant flux 
are all above 32%. 

We conclude that the VHE fluxes are constant on timescales 
of 30 — 75 minutes within the MAGIC sensitivity. 

4.2. Search for periodicity 

The emission of the binary system changes periodically in 
radio, optical and X-rays, and the modulation is associated 



with the orbital period of the binary system. At higher ener- 
gies, EG RET measurements showed hints for variable 7-ray 
emission (Tava nTet al.lll998l) . although no periodicity could 
be established with these data. One of the aims of the MAGIC 
long observational campaign on LS I +61°303 was to search 
for periodic VHE 7-ray emission. 

In order to maximize the sensitivity and accuracy of the tim- 
ing analysis we used the data presented i n this work togethe r 
with the data taken in the first campaign (AlberteLalJ|2006), 
with observation time of 54 h and covering 6 orbital cycles. 

The periodicity analysis was carried out using the formal- 
ism d eveloped by Lomb and Scargle ( Lombl [T 976t fScargle 
1982). This formalism allows to analyze unevenly sampled 
data while still keeping the simple exponential probability 
density function (PDF, P(z > zq) — e~ z °) for Gaussian 
White Noise (GWN) as valid for the classical Fourier analysis 
of evenly sampled data. A remaining problem caused by the 
uneven data sampling is that the independent Fourier spacing 
is broken, i.e. even a single frequency component can result 
in a complex power spectrum with a large number of aliasing 
peaks. An other problem is that the mean and variance values 
that enter the Lomb-Scargle periodogram have to be estimated 
from the data themselves. 

A practical method to determine the chan ce probabilities is 
the following (see e.g. lFrescura et alj (120071) ): 

1. A large number of random data series is constructed 
with a Monte Carlo simulation of random fluxes while 
keeping the sampling times fixed. 

2. For each random series, we construct a periodogram, 
sampling it for a pre-selected group of frequencies. 

3. For each frequency, we compare the periodogram de- 
rived from the real data set with the probability density 
function (PDF) obtained from the simulated random se- 
ries, in order to empirically determine the (pre-trial) 
chance probability. 

4. The overall (post-trial) chance probability is computed 
according to the following generalization: for each sim- 
ulated data series we inspect the corresponding peri- 
odogram, identify the highest Fourier power that occurs 
at any of the pre-selected frequencies, and use this value 
to construct the post-trial PDF. It should be noted that 
this constructed PDF is based on the null hypothesis of 
GWN. 

Integration of the post-trial PDF gives the complementary Cu- 
mulative Distribution Function (cCDF) which is used to de- 
termine the (post-trial) chance probability for a given Fourier 
power value. 

In Fig. [5] we show the empirical post-trial cCDF of the 
Lomb-Scargle power, estimated via Monte Carlo simulation 
of random fluxes. The expected cCDF above a spectral peak 
z Q is F(z > z ) = 1 - (1 - e~ z °) M , where M is the 
number of independent frequencies. By fitting the PDF for 
LS I +61°303, we obtain a probability of 75% (x 2 /dof = 
263.9/279) and a number of independent frequencies of M = 
550.8 ± 0.6. This result is used to estimate the chance proba- 
bility of the Lomb-Scargle powers. 

In Fig. [6] (middle panel) we show the Lomb-Scargle pe- 
riodogram for an almost (up to detector related effects) in- 
dependent background sample, obtained simultaneously with 
the LS I +61°303 data (see below for the time intervals). The 
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FIG. 5. — Post-ti'ial complementary Cumulative Distribution Function for 
the Lomb-Scargle power derived from 10 e random time series. The expected 
cCDF is also indicated (solid line). This function is used to estimate the 
chance probability for powers above 20. 
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FIG. 6. — Lomb-Scargle periodogram over the combined 2005 and 2006 
campaigns of LS I +6 1 ° 303 data (upper panel) and simultaneous background 
data (middle panel). In the lower panel we show the periodograms after sub- 
traction of a sinusoidal signal (see Fig. [7J at the orbital period (yellow line) 
and a sinusoidal plus a Gaussian wave form (blue line). Vertical dashed line 
corresponds to the orbital frequency. Inset: zoom around the highest peak, 
which corresponds to the orbital frequency (0.0377d — 1 ). Its post-trial prob- 
ability is nearly 10 -7 (see Fig. [5}. The IFS is also shown. 



highest obtained power is 7.5, which yields a probability of 
0.3. Thus we obtain no significant probability peaks for any 
of the scanned frequencies. 

We apply the Lomb-Scargle test to the LS I +61°303 data 
and obtain the periodogram shown in Fig. [6] (upper panel). 
The periodogram is performed with the LS I +6 1 303 integral 
flux above 400 GeV, measured in a time interval [ij — =r , tj + 

for At = 15 minutes and i = 0, . . . 717 data points. The 
overall time range is 442 days, which yields an independent 
Fourier spacing (IFS) of v IFS = 1/T = 0.0023 d" 1 . We 
scanned the frequency range 0.0023-0.25 d _1 with a an over- 
sampling factor of 5. 
A maximum peak in the Lomb-Scargle periodogram is 




0.8 1 
Orbital phase 

FIG. 7.— LS I +61°303 7-ray flux above 400 GeV obtained from the 
first and second campaigns, folded with the orbital frequency in bins of 0.05 
in phase. The black curve is a fit to a sinusoidal signal. We also fitted a 
sinusoidal signal plus a Gaussian component (blue dotted line), which ad- 
justs better to the data (fit parameters are given in the inset). The vertical 
dashed lines mark the two me asurements of the p eriastron passage (accord- 
ing to Casares et al. 1 2005) and Grundstrom et al. l 2003)). 



clearly seen at frequency v = 0.0373d -1 , for which we ob- 
tain a power of ~ 22, corresponding to a post-trial chance 
probability of 2 x 10~ 7 . 

Several less prominent but significant peaks are also de- 
tected for other frequencies (e.g. 0.041 d _1 with probability 
< 10~ 5 ). Those peaks are related to the signal, since they 
are not present in the contemporaneous background sample 
(Fig- IU middle panel). These are aliasing peaks of the orbital 
period of LS I +61°303 caused by the various gaps in the data 
set. 

The observational bias due to the moon cycle cannot be the 
responsible for the observed peak since this period should oth- 
erwise be also present in the background periodogram. 

The data folded with the peak frequency (v — 0.0377d _1 ) 
is presented in Fig. [7] where a sinusoidal fit is performed 
(x 2 /dof = 123.6/16). Subtracting the obtained sinusoid 
from the data, we produce the periodogram shown in Fig. [6] 
(lower panel, yellow line). The peak associated to the orbital 
frequency has been removed as expected. Also the satellite 
peaks are reduced, but the fact that some of the other peaks 
do not achieve a level consistent with the background test in- 
dicates that the signal in the LS I +61°303 data is not purely 
sinusoidal. 

To reduce these remaining powers, we fitted the data set 
with a more complex signal. Motivated by the data shape, we 
fitted the data set with a sinusoidal function plus a Gaussian 
signal contribution (x 2 /dof = 58.1/13), as shown in Fig. [7] 
(blue dotted line). The corresponding periodogram subtract- 
ing this function to the LS I +61°303 data set is given in Fig. [6] 
(lower panel, blue line). The orbital frequency peak has been 
removed and some of the periodogram peaks are much more 
reduced than in the purely sinusoidal subtraction, giving a bet- 
ter agreement with the background periodogram level. 

We performed a Monte Carlo simulation to evaluate the er- 
ror in the frequency estimation without any signal shape as- 
sumption: we simulate light curves where the number of 7- 
ray and background candidates are selected randomly from 
Poisson distributions with a mean equal to the actually mea- 
sured distributions of events, arriving in every given time in- 
terval. The periodogram is calculated for 10 3 of those ran- 
domly generated series, and the distribution of the result- 
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FIG. 8. — LS I +61°303 period measurements in different wavelengths. 
Blue band indicates a 3a region around the radio measurement. The 7-ray 
period (in red) is compatible within 1.5<r with it. 

ing peak power frequencies is fitted with a Gaussian func- 
tion, yielding an error of 0.0003 d _1 . An accurate peak fre- 
quency determination is done by scanning more frequencies 
(increasing the oversampling factor) around the frequency 
which has maximum probability in the periodogram, and is 
found to be (0.0373±0.0003) d~\ corresponding to a period 
of (26.8±0.2)days. 

In Fig. |8] we show the period obtained with MAGIC 
data compared to the measurements in other wavelengths. 
The most accurate measure of the or bital period i s 
(26.4960±0.0028) days, reported in radio bv lGregory|(l2002l) . 
We also show period measuremen t s repo rted in near IR 
and opt ical V-band dParedes et al.l 1 19941) . optical wave- 
lengths dMendelson and Mazehl 1 19891). photo metry in the I- 
B and I bands ( Mendelson an d Mazehl 1 19941) . Ha measure- 
ments (Zamano v et alJ[l999|), and soft X-ray m easurements 
from lWen et al I (l2006h and lParedes et alJ (1 19971) . The period 
obtained with MAGIC data is compatible with the orbital ra- 
dio measurement within 1.5er. 

5. CONCLUSION 

We find that LS I +61°303 is a periodic 7-ray binary with 
an orbital period of 26.8±0.2 days (and chance probability 



~ 10 ), compatible with the optical, radio and X-ray period. 
This result implies that the flux modulation is tied to the or- 
bital period. The high state in VHE 7-rays occurs in the same 
phases as the X-ray high state. This is especially interesting 
since we found a additional hint for X-ray/7-ray variability 
correlation in the orbital phase 0.85. A strictly simultaneous 
multi-wavelength campaign is needed to investigate this cor- 
relation in more detail. 

We looked for possible intranight variability and found the 
flux consistent with being constant within errors in 30-75 min- 
utes time-scales. 

We produce energy spectra for two phase bins 0.5 < (j> < 
0.6 and 0.6 < cf) < 0.7 and averaged flux values for several 
phase bins. There is clear evidence for a significant change 
in the VHE 7-ray flux level between different phase bins of 
LS I +61°303. The spectral photon index does not show this 
dependence on the phase. All derived spectral photon indices 
are compatible with 2.6 ± 0.2 , obtained from the most signif- 
icant phase bin of LS I +6 1 303 . 

We can put constraints to the emission at the periastron pas- 
sage and conclude that the system is detected in 7-rays only 
in the phases 0.5 — 0.9. Since significant emission is only 
detected in an orbital sector off the phases at which the max- 
imum gamma ray flux should occur under photon-photon ab- 
sorption (see fig. 5 in iDubusI 12006), the latter can hardly be 
the only source of variability in the emission. 

Thorough multiwavelength observations will allow us to 
probe the intrinsic variability of the non-thermal emission 
from LS I +61°303 along the orbit and can proof possible 
correlations between the X-ray and TeV energy bands. This 
is a necessary step for understanding the source nature, and 
the physics underlying the VHE radiation. 
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